# https://blog.csdn.net/XiongSiqi_blog/article/details/130978226
import libnum
import sys

def HGCD(a, b):
    if 2 * b.degree() <= a.degree() or a.degree() == 1:
        return 1, 0, 0, 1
    m = a.degree() // 2
    a_top, a_bot = a.quo_rem(x^m)
    b_top, b_bot = b.quo_rem(x^m)
    R00, R01, R10, R11 = HGCD(a_top, b_top)
    c = R00 * a + R01 * b
    d = R10 * a + R11 * b
    q, e = c.quo_rem(d)
    d_top, d_bot = d.quo_rem(x^(m // 2))
    e_top, e_bot = e.quo_rem(x^(m // 2))
    S00, S01, S10, S11 = HGCD(d_top, e_top)
    RET00 = S01 * R00 + (S00 - q * S01) * R10
    RET01 = S01 * R01 + (S00 - q * S01) * R11
    RET10 = S11 * R00 + (S10 - q * S11) * R10
    RET11 = S11 * R01 + (S10 - q * S11) * R11
    return RET00, RET01, RET10, RET11

def GCD(a, b):
    print(a.degree(), b.degree())
    q, r = a.quo_rem(b)
    if r == 0:
        return b
    R00, R01, R10, R11 = HGCD(a, b)
    c = R00 * a + R01 * b
    d = R10 * a + R11 * b
    if d == 0:
        return c.monic()
    q, r = c.quo_rem(d)
    if r == 0:
        return d
    return GCD(d, r)

sys.setrecursionlimit(500000)

n = 128134155200900363557361770121648236747559663738591418041443861545561451885335858854359771414605640612993903005548718875328893717909535447866152704351924465716196738696788273375424835753379386427253243854791810104120869379525507986270383750499650286106684249027984675067236382543612917882024145261815608895379
e = 5331
c1 = 60668946079423190709851484247433853783238381043211713258950336572392573192737047470465310272448083514859509629066647300714425946282732774440406261265802652068183263460022257056016974572472905555413226634497579807277440653563498768557112618320828785438180460624890479311538368514262550081582173264168580537990
c2 = 43064371535146610786202813736674368618250034274768737857627872777051745883780468417199551751374395264039179171708712686651485125338422911633961121202567788447108712022481564453759980969777219700870458940189456782517037780321026907310930696608923940135664565796997158295530735831680955376342697203313901005151

a=2
b=1

R.<x> = PolynomialRing(Zmod(n))
f = x^e - c2
g = (a*x+b)^e - c1

res = GCD(f,g)

tmp=((int(-res.monic().coefficients()[0])))
print(tmp)
m=tmp*2+1
import libnum
print(libnum.n2s(int(m)))